AD-A1S0  036 


File  cup; 


REPORT  DOCUMENTATION  PAGE 

lib.  RESTRICTIVE  MARKINGS 


2b.  DECLASSIFICATION  /  DOWNUKA 


4.  PERFORMING  ORGANIZATION  RE? 


6«.  NAME  OF  PERFORMING  ORGANIZATION 
Department  of  Systems  Science 
and  Mathematics 


6c  ADDRESS  (Gty,  State,  and  ZIP  Code) 

Washington  University 

Campus  Box  1040 
St.  Louis,  MO  63130 


8a.  NAME  OF  FUNDING  /SPONSORING 
ORGANIZATION  AFOSR/NM 


8c  ADDRESS  (City,  State,  and  ZIP  Code) 

Building  410 

Bolling  AFB ,  DC  20332-6448 


11.  TITLE  (Include  Security  Classification) 


3.  DISTRIBUTION /AVAILABILITY  OF  REPORT 
Distribution  of  this  report  is  unlimited 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

AFOSR-TR.  8  8-0 


7a.  NAME  OF  MONITORING  ORGANIZATION 
AFOSR/PKZA 


7b  ADORESS  (City,  State,  and  ZIP  Code) 

Building  410 

Bolling  Air  Force  Base,  DC  20332-6448 


8b  OFFICE  SYM80L  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

*  AFOSR-82-0315  C.  D 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

NO 

2304 


WORK  UNIT 
ACCESSION  NO 


DEVELOPMENT  AND  APPLICATION  OF  THE  P-VERSION  OF  THE  FINITE  ELEMENT  METHOD 


12  PERSONAL  AUTHOR(S) 


13a.  TYPE  OF  REPORT 


I.  Norman  Katz 


13b.  TIME  COVERED 


Final  ScientiftcfROM 


30  to  87/9/30 


14.  DATE  OF  REPORT  (Year,  Month,  Day)  |l  5  PAGE  COUNT 

87/12/30  27 


COSATI  COOES 


SUB-GROUP 


18.  SUBJECT  TERMS  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 


19.  ABSTRACT  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

' ^  The  p-version  of  the  finite  element  method  is  a  new,  important,  computationally 
efficient,  approach  to  finite  element  analysis.  It  is  more  robust  than  the  conventional 
h-version  and  its  rate  of  convergence,  for  domains  with  corners  and  for  other  singularity 
problems,  is  twice  that  of  the  h-version. 

Hierarchic  elements  which  implement  the  p-version  efficiently  have  been  formulated 

so  as  to  enforce  C  or  C  continuity  in  the  planar  case,  and  so  as  to  enforce  C^  con¬ 
tinuity  in  three  dimensions. 

*  Continued  on  the  reverse  side 


20  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

OUNCLASSIFIEQ/UNLIMITED  □  SAME  AS  RPT  Q  OTIC  USERS  UNCLASSIFIED 


22c  OFFICE  SYMBOL 

m 


m 


DO  FORM  1473, 84  mar 


83  APR  edition  may  be  used  until  exhausted. 
All  other  editions  are  obsolete. 


SECURITY  CLASSIFICATION  OF  ~HIS  PAGE 
r  * 

!  j> 

....  it  I 


Recent  research  accomplishments  include: 

1.  Development  of  an  algorithm  that  finds  all  roots  of  an  analytic  func¬ 
tion  in  a  finite  domain. 

2.  Preprocessing  procedures  to  restrict  the  search  in  unbounded  domains 
which  contain  roots  to  bounded  domains, 

3.  A  reliable  numerical  argument  principle  algorithm  to  compute  number  of 
zeros  within  a  closed  contour. 

4.  Formulation  of  equations  which  determine  the  nature  of  stress  sin¬ 
gularity  at  a  corner  of  a  plate  composed  of  n  isotripic  materials. 

All  of  the  above  are  used  in  the  extraction  method  for  p-version  finite 
element  analysis  of  composite  materials  with  corners. 


V  i 


AiO^K-TK-  b8-0l4y 


Annual  and  Final  Scientific  Report 

Air  Force  Office  of  Scientific  Research  Grant  AFOSR  82-0315  C,  D 

Period:  30  September  1985  -  30  September  1987 

Title  of  Research:  Development  and  Application  of  the 

p-version  of  the  Finite  Element  Method 

Principal  Investigator:  I.  Norman  Katz,  Professor  of 

Applied  Mathematics  and  Systems 
Science 


Co-principal  Investigator:  Barna  A.  Szabo,  A.P.  Greensfelder 
(  Professor  of  Mechanical  Engineering 


<1 


Department  of  Systems  Science  and  Mathematics 
School  of  Engineering  and  Applied  Science 
Washington  University  -  Campus  Box  1040 
St.  Louis,  Missouri  63130 


vh*  v;v:v:ww*v. 


TABLE  OF  CONTENTS 


Introduction . ' . 

Summary  of  Earlier  Research  Accomplishments  .  .  . 

2.1  Rate  of  Convergence  of  the  p -version  .  .  .  , 

2.2  Hierarchic  Families  of  Solid  Finite  Elements 


2.3  Quantities  of  Special  Interest 

2.3.1  Use  of  Functional  Forms 

2.3.2  Extraction  Techniques 


Summary  of  Recent  Research  Accomplishments 


3.1  Extraction  Techniques  for  Composite  Materials  - 
Formulation  of  Equations  . 


3.2  Globally  Convergent  Algorithms  for 
Finding  all  the  Eigenvalues  .  .  .  . 


3.3  A  Numerical  Argument  Principle  . 

Computer  Implementations  and  Technology  Transfer  . 

4.1  COMET-X  . 

4.2  FIESTA/3D  . 

4 . 3  PROBE  . 

Figures  . 

References  . 

Principal  Personnel  . 

Papers  Published  and  Presented  Since  Start  of  the  Project  .  .  . 

8.1  Published  Papers  . 

8.2  Presented  Papers  . 

8.3  Visits  and  Seminars  at  Government  Laboratories  . 


'  <\ 

w  * 

<Xy 


For 
26  I 


j  li-  i  C  7  A.  ri 
1  Unannounced 
!  Just  If  i'"U  i ->n . 


By  .  .  _ 

Di  3t  i-ibut  l  >»/ 

Ava liability  Codes 
Avail  tuui/or 


;oist 

Sptc 

|Hj 

1. 


INTRODUCTION 


There  are  now  three  basic  approaches  to  finite  element  analysis.  In 
all  approaches  the  domain  fi  is  divided  into  simple  convex  subdomains 
(usually  triangles  or  rectangles  in  two  dimensions,  and  tetrahedra  or  bricks 
in  three  dimensions)  and  over  each  subdomain  the  unknown  is  approximated  by 
a  (local)  basis  function  (usually  a  polynomial  of  degree  <  p) .  Basis  func¬ 
tions  are  required  to  meet  continuously  at  boundaries  of  subdomains  in  the 
case  of  planar  or  3  dimensional  elasticity,  or  smoothly  in  the  case  of 
plate  bending.  The  approaches  are: 

1.  The  h-version  of  the  finite  element  method.  In  this  approach  the 
degree  p  of  the  approximating  polynomial  is  kept  fixed,  usually  at 
some  low  number  such  as  2  or  3.  Convergence  is  achieved  by  allowing 
h,  the  maximum  diameter  of  the  convex  subdomains,  to  go  to  zero. 
Estimates  for  the  error  in  energy  have  long  been  known  [1,  2,  3].  In 
all  of  these  estimates  p  is  assumed  to  be  fixed  and  the  error  es¬ 
timate  is  asymptotic  in  h,  as  h  goes  to  zero. 

2.  The  p-version  of  the  finite  element  method.  In  this  approach  the 
subdivision  of  the  domain  Q  is  kept  fixed  but  p  is  allowed  to 
increase  until  a  desired  accuracy  is  attained.  The  p-version  is 
reminiscent  of  the  Ritz  method  for  solving  partial  differential  equa¬ 
tions  but  with  a  crucial  distinction  between  the  two  methods.  In  the 
Ritz  method  a  single  polynomial  approximation  is  used  over  the  entire 
domain  (n,  in  general,  is  not  convex).  In  the  p-version  of  the 
finite  element  method  polynomials  are  used  as  approximations  over 
convex  subdomains.  This  critical  difference  gives  the  p-version  a  more 
rapid  rate  of  convergence  than  either  the  Ritz  method  or  the  h-version. 
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3.  The  h-p  version  of  the  finite  element  method.  In  this  approach  both 

the  degree  p  of  the  approximating  polynomial  and  the  maximum  diameter 
h  of  the  convex  subdomains  are  allowed  to  change. 

The  p-version  of  the  finite  element  method  requires  families  of  polyno¬ 
mials  of  arbitrary  degree  p  defined  over  different  geometric  shapes. 
Polynomials  defined  over  neighboring  elements  join  either  continuously  (are 

in  C^)  for  planar  or  three  dimensional  elasticity,  and  smoothly  (are  in 

C^")  for  plate  bending.  In  order  to  implement  the  p-version  efficiently  on 
the  computer,  these  families  should  have  the  property  that  computations 
performed  for  an  approximation  of  degree  p  are  re-usable  for  computations 
performed  for  the  next  approximation  of  degree  p  +  1.  We  call  families 
possessing  this  property  hierarchic  families  of  finite  elements. 

The  h-version  of  the  finite  element  method  has  been  the  subject  of 
intensive  study  since  the  early  1950' s  and  perhaps  even  earlier.  Study  of 
the  p-version  of  the  finite  element  method,  on  the  other  hand,  began  at 
Washington  University  in  St.  Louis  in  the  early  1970 's  and  led  to  a  more 
recent  study  of  the  h-p  version.  Research  in  the  p-version  (formerly  called 
The  Constraint  Method)  has  been  supported  in  part  of  the  Air  Force  Office  of 
Scientific  Research  since  1976. 


2.  SUMMARY  OF  EARLIER  RESEARCH  ACCOMPLISHMENTS 
2 . 1  Rate  of  Convergence  of  the  p-version 

Extensive  computational  experiments  have  long  furnished  empirical 
evidence  that  the  rate  of  convergence  of  the  p-version  is  significantly 
higher  than  that  of  the  h-version  (see  for  example  [4,  5,  6]). 

The  two  theorems  described  in  detail  in  [5,  6]  provide  a  mathematical 
explanation  for  the  efficiency  of  the  p-version. 


Roughly  speaking,  these  theorems  state  that  if  the  domain  is  a  bounded 
polygon  in  two  dimensions,  and  if  the  criterion  used  to  measure  efficiency 
is  the  number  of  degrees  of  freedom  required  to  achieve  a  given  error  in 
energy,  then  the  rate  of  convergency  of  the  p-version  is  twice  that  of  the 
h-version. 

Also,  it  is  shown  in  [7,  8]  if  optimal  mesh  refinement  (not  necessarily 
quasi -uniform)  is  combined  with  optimal  p-distribution  then  the  asymptotic 
convergence  is  exponential. 

Figure  1  illustrates  graphically  the  rates  of  convergence  of  the  h- , 
p- ,  and  h-p-versions . 

2 . 2  Hierarchic  Families  of  Solid  Finite  Elements 

In  order  to  implement  the  p-version  efficiently,  families  of  finite 
elements  are  needed  with  the  hierarchic  property:  computations  performed 
for  an  approximation  of  order  p  should  be  re-usable  when  raising  the  order 
to  p+1.  More  specifically,  the  stiffness  matrix  corresponding  to  the 
polynomial  approximation  of  degree  p  should  be  a  submatrix  of  the  polyno¬ 
mial  approximation  of  degree  p+1.  In  terms  of  basis  functions,  this 
implies  that  the  basis  functions  for  a  pth  order  approximation  should  be  a 
subset  of  the  basis  functions  for  a  (p  +  l)st  order  approximation. 

Hierarchic  families  for  triangles  both  in  the  and  in  the  case 

are  described  in  detail  in  [9,  10,  11,  12,  13,  14], 

2 . 3  Quantities  of  Special  Interest  -  Extraction  Techniques 

Often  the  main  purpose  of  a  finite  element  analysis  is  to  obtain  values 


of  a  few  important  quantities  with  a  high  accuracy.  In  structural 
mechanics,  for  example,  the  values  of  displacements  or  stresses  in  a  small 
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number  of  designated  areas,  or  the  stress  intensity  factor  at  a  small  number 
of  points  is  of  critical  importance  for  design.  Both  the  h-  and  p- 
versions  of  the  finite  element  method  give  approximations  for  these  values. 
However,  it  is  much  more  efficient  to  use  a  post-processing  technique  which 
uses  weighted  averages  of  values  taken  directly  from  the  finite  element 
approximation.  The  post-processor  determines  these  quantities  of  special 
interest  much  more  accurately.  In  particular,  when  stresses  are  computed 
pointwise  as  derivatives  of  displacements  in  the  p-version,  they  may  exhibit 
(some  times  severe)  oscillatory  behavior.  In  the  case  of  the  centrally 
cracked  panel  shown  in  Figure  2(h),  for  example  the  normal  stresses  a 

along  the  x-axis  computed  for  polynomial  orders  ranging  from  1  to  7  are 

shown.  The  oscillatory  behavior  near  the  crack  tip  singularity  is  evident. 
Two  techniques  are  being  developed  for  the  post-processing  of  quantities  of 
special  interest.  Both  techniques  are  very  accurate,  yielding  approxima¬ 
tions  that  are  of  the  same  order  of  accuracy  as  the  strain  energy.  (This  is 

the  square  of  the  error  in  energy  norm) . 

2.3.1  Use  of  Functional  Forms 

This  approach  is  based  on  the  idea  that  the  functional  forms  of  the 

quantities  of  interest  are  generally  known.  In  the  case  of  a  centrally 

cracked  panel,  for  example,  the  displacements  u(z) ,  v(z) ,  and  the  stresses 

a  (z),  t  (z)  and  a  (z)  (z  -  x  +  iy)  are  given  in  terms  of  two  functions, 
x  xy  y 
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j-0  J 


are  holomorphic  functions.  Approximate  values  for  the  coefficients  a  , 
are  determined  using  the  displacements  u^(z),  v^(z)  computed  by  the  p- 

version  of  the  finite  element  method. 

This  technique  was  used  in  [15]  to  obtain  improved  values  of  a ^  along 

the  x-axis  in  the  case  of  a  centrally  cracked  panel.  Sample  results  are 
shown  in  Figure  3  (R  in  Figure  3  is  the  radius  of  the  circle  surrounding  the 
crack  tip  within  which  values  of  up(z)  and  vp(z)  were  taken) .  Figure  4 

shows  the  convergence  of  this  post  processing  technique  for  the  stress 


intensity  factor  K  as  a  function  of  1/NDF. 


2.3.2  Extraction  Techniques 

A  general  form  for  post  processing  calculations  is  given  in  [16,  17, 
18].  Green's  theorem  and  generalized  influence  functions  are  used  together 
with  smooth  cut-off  functions  and  blending  functions  in  order  to  calculate 
higher  derivatives  of  the  unknown  function  and  also  stress  intensity 
factors . 


i 


I 

I 
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These  techniques  have  been  applied  to  post  processing  of  the  u  and 


Vp  displacement  fields  obtained  from  the  p-version. 
[19,  20], 


See,  for  example, 


3.  SUMMARY  OF  RECENT  RESEARCH  ACCOMPLISHMENTS 

3 . 1  Extraction  Techniques  for  Composite  Materials-Formulation  of  Equations 
In  the  techniques  discussed  in  2.3.1  and  2.3.2  for  the  post-process 
calculation  of  stresses  and  stress  intensity  factors  it  is  essential  that 
the  behavior  of  the  solution  at  the  singular  point  be  known  in  advance.  It 
is  well  known  [21]  that  the  behavior  of  the  displacements  of  a  thin  plate  in 
a  state  of  plain  stress  or  plain  strain  at  a  vertex  of  a  sector  is  of  the 
form 


I  VAk  los  q* 

where  (r ,0)  are  polar  coordinates  centered  at  the  vertex  and  0k(0)  are 

smooth  functions  (usually  sines  or  cosines).  The  (complex)  numbers  A^,  and 

their  multiplicities  are  determined  by  enforcing  boundary  conditions  at  both 
sides  meeting  at  the  vertex.  The  values  of  A^  with  smallest  real  part 


have  been  determined  in  a  well  known  paper  by  Williams  [22]  for  different 
sets  of  boundary  conditions  in  the  case  that  the  sector  is  composed  of  a 
single  material. 

In  many  modern  applications  the  sector  is  composed  of  N  materials 
each  with  its  own  Poisson's  ratio  v ^  and  Young's  modulus  E^.  In  this 

case  the  determination  of  the  eigenvalues  A^  and  corresponding  eigenfunc¬ 


tions  r^k  log  ^k  0^(0) 


is  considerably  more  difficult.  However,  it  is 


V.' 
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important  to  determine  the  eigenvalues  A^  and  corresponding  eigenfunctions 

in  order  to  use  the  post  processing  capabilities  of  the  p-version. 

Consider  N  materials  meeting  at  a  point  as  shown  in  the  adjoining 
figure.  The  interface  conditions  at  edges  1  through  N-l  are:  edge  N 


rr9,i(*i)  ~  rr 9,  i+1  (fli) 


a9,  i  (V  "  °9,  i+1  (*i> 


Yi  “Y  i+i  (*i> 


U  .  (9.)  -  U  .  (9.) 

r,i  v  i'  r,  i+i  v  i' 


Let  n ^  — 


the  shear  modulus  in  material  i 


2(l+v1) 
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Hi  " 


be  the  unknown  e 
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and  define  the  boundary  matrices  Bj  for  j  -  0,  N  by 

Bq  “  T2  (free  edge)  B^  -  T2  (free  edge) 


Bq  -  T2  +  F  T3  (clamped  edge) 


BN  "  T2  +  F  T3  (clamped  edge) 


It  was  shown  [23,  24]  that  each  eigenvalue  A  satisfies  det  S  -  0 
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Furthermore  a  lengthy  computation  shows  that 


N  N- 1  -  2 

det  S  -  det  S"  (A+l)  (-64)1  L  n  EV 

t-2 
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Note  that  this  last  expression  gives  S"  as  product  of  a  2x4  matrix, 
followed  by  (4x4)  matrices,  followed  by  a  (4x2)  matrix.  Therefore  S"  is  a 
(2x2)  matrix. 

It  is  of  course  sufficient  to  solve  det  S"  -  0.  The  expression  for 
S"  even  in  the  simplified  form  given  above  is  vary  lengthy.  The  symbolic 
manipulator  code  REDUCE  can  be  used  to  obtain  the  expression  for  S".  A 
closed  from  would  be  much  more  efficient  for  numerical  computations  and  is 
currently  being  sought. 

E. 

It  is  remarkable  that  if  the  shear  modulus  /i.  -  -  -  n  is  con- 

1  2(1+*.) 

stant  for  all  materials  then  a  relatively  simple  closed  form  nonlinear 
equation  can  be  obtained  for  each  eigenvalue  A: 


ft  N 
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Once  the  nonlinear  equations  for  the  eigenvalues  A  are  formulated, 
they  are  solved  for  eigenvalues  with  smallest  real  parts.  These  eigen¬ 
values,  together  with  associated  eigenfunctions  are  used  in  the  extraction 
techniques  discussed  earlier  in  order  to  determine  the  constants  a^.  The 

displacement  and  stress  fields  in  the  neighborhood  of  a  vertex  formed  by 
composite  materials  are  then  computable  through  the  post  processor  of  the  p- 
version.  This  provides  a  powerful  new  procedure  for  the  possible  prediction 
of  crack  initiation  based  on  factors  (the  a^)  which  are  analogous  to  stress 

intensity  factors  (used  for  the  prediction  of  crack  propagation) . 

3 . 2  Globally  Convergent  Algorithms  for  Finding  all  the  Eigenvalues  X 

The  determination  of  all  of  the  eigenvalues  A^  themselves  and  their 

corresponding  eigenfunctions  is  a  challenging  problem.  We  have  developed 
algorithms  with  the  following  properties: 

1.  All  eigenvalues  in  a  given  infinite  strip  -Xq  <  ReA  <  Xq  are 

computed  by  the  algorithms  to  a  prespecified  accuracy. 

2.  All  of  the  corresponding  eigenfunctions  are  computed. 

It  is  important  that  no  eigenvalues  (with  small  real  parts)  be  omitted 
in  the  solution  of  the  eigenequation  because  this  will  lead  to  incorrect 
characterization  of  the  singularity  at  the  corner.  The  extraction  technique 
described  earlier  will  then  fail. 

For  this  purpose  we  have  developed  globally  convergent  algorithms. 

Our  approach  is  based  on  an  exclusion  method.  In  this  approach,  a  simple 
criterion  is  formulated  to  determine  whether  a  root  does  not  lie  in  a  given 
region.  If  so,  the  region  is  excluded  from  further  computation.  This 
approach  assures  that  all  roots  of  a  nonlinear  function  in  a  given  starting 


region  will  be  determined.  Then,  those  with  smallest  real  part,  can  be 
selected. 


In  [25],  a  reliable  solver  for  all  roots  of  a  nonlinear  equation  is 
described  in  detail.  The  idea  of  using  global  information  (in  the  form  of 
bounds  on  derivatives)  dynamically  is  introduced.  Also  the  concepts  of 
dominating  function  and  dominated  functions  are  formulated  to  reduce  the 
search  domain  from  an  infinite  strip  to  a  bounded  domain. 

3 . 3  A  Numerical  Argument  Principle 

In  order  to  assure  that  all  roots  have  indeed  been  obtained,  we  have 
studied  the  use  of  the  principle  of  the  argument  for  functions  of  a  complex 
variable,  to  determine  the  number  of  zeros  of  the  eigenequations  in  a  given 
(bounded)  domain.  It  is  crucial  that  the  steps  taken  in  numerically  in¬ 
tegrating  along  a  contour  be  small  enough  to  assure  that  there  are  no 
spurious  jumps  in  the  change  of  the  argument  which  may  be  caused  by 
inaccuracies.  A  criterion  has  been  developed  to  compute  the  size  of  these 
steps.  Details  are  given  in  [26], 

These  algorithms  may  have  wider  application  to  the  problem  of  finding 
all  zeros  of  anv  analytic  functions  in  a  given  domain. 

4.  COMPUTER  IMPLEMENTATIONS  AND  TECHNOLOGY  TRANSFER 

Our  earlier  work  on  research  and  development  of  the  p-version  have 
resulted  in  several  computer  implementations.  As  our  research  progresses, 
new  procedures  and  ideas  are  continually  incorporated  into  these  codes 
making  them  more  efficient.  Two  of  the  codes  (FIESTA/3D  and  PROBE)  are  in 
commercial  use.  A  detailed  description  of  the  efficient  technology  transfe 


was  reported  in  SIAM  NEWS  Vol.  19,  No.  6,  (November  1986)  and  Vol.  20,  No.  2 
(March  1987) . 

4.1  COMET -X 

COMET-X  is  an  experimental  computer  code  which  implements  the  p -vers ion 
of  the  finite  element  method  by  using  the  hierarchic  families  which  have 
been  constructed.  COMET-X  is  maintained  by  the  Center  for  Computational 
Mechanics  at  Washington  University.  COMET-X  can  be  used  as  a  code  to  imple¬ 
ment  the  h-version  as  well  simply  by  fixing  the  polynomial  order  p  and 
refining  the  mesh. 

COMET-X  currently  has  the  following  capabilities: 

A.  Element  types:  Stiffeners,  triangular  elements,  triangular  elements 
with  one  side  curved,  rectangular  elements,  solid  elements  of  the  shapes 
described  earlier. 

B.  Types  of  Analysis:  Laplace  and  Poisson  equations,  plane  elasticity, 
temperature  distribution  in  3-dimensions. 

C.  Special  Capabilities:  non-uniform  p-distribution,  elastic  fracture 
mechanics  computations  in  two  dimensions,  nearly  incompressible  solids, 
linear  boundary  layer  problems. 

D.  Pre-  and  Post  processing  capabilities  including  graphics  and  visual 
displays.  The  capabilities  and  usage  of  COMET-X  are  described  in  detail 
in  [27]. 


4.2  FIESTA/3D 

FIESTA/3D  is  a  software  system  for  static  analysis  of  solid  structures 
based  on  the  p-version  of  the  finite  element  method.  It  is  marketed  by 
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McDonnell -Douglas  Automation  Company.  Some  of  the  advanced  features  avail¬ 
able  on  FIESTA/3D  are: 

•  controllability  of  the  quality  of  the  solution 

•  efficient  modeling 

•  treatment  of  stress  singularities 

•  a  method  for  surface  identification. 

4.3  PROBE 

PROBE  is  an  advanced  computer  implementation  of  the  p -vers ion  for  2- 
dimensional  analysis  currently  under  development  at  NOETIC  Technologies  in 
St.  Louis,  Missouri.  Some  of  PROBE'S  unique  features  are: 

•  automatic  error  estimators.  These  estimators  provide  initial  feed¬ 
back  on  solution  quality  by  computing  element-by-element  and  edge- 
by-edge  equilibrium  checks. 

•  interactive  error  estimators.  These  estimates  offer  immediate 
feedback  on  solution  quality  by  computing  action/reaction  checks 
overall  equilibrium  checks  and  a  convergence  trajectory  of  stresses 
and  strains  at  user-selected  points. 

•  load  tracking,  for  extracting  free -body  diagrams 

•  pinpoint  solutions  which  give  specific  results  anywhere  in  the  model 

•  Mode  I  (Symmetric)  and  Mode  II  (antisymmetric)  Stress  Intensity 
Factors  for  Fracture  Mechanics 

•  Precise  curve  definitions  which  eliminate  mesh  refinement  near 
cutouts  and  provide  more  precise  solutions  in  critical  regions. 


6. 
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July  16-20,  1984. 

"Implementation  of  a  Finite  Element  Software  System  with  h-  and  p- 
extension  Capabilities,"  by  B.A.  Szabo,  8th  Invitational  Symposium  on 
the  Unification  of  Finite  Element- Finite  Differences  and  the  Calculus 
of  Variations,  University  of  Connecticut,  Storrs,  Connecticut,  May  3, 
1985. 

"Stress  Singularities  at  Angular  Corners  of  Composite  Plates,"  X-R. 

Ying  and  I.N.  Katz,  SIAM  Spring  Meeting,  Pittsburgh,  PA,  June  1985. 

"Computation  of  the  Amplitude  of  Stress  Singular  Terms  for  Cracks,  and 
Reentrant  Corners,"  B.A.  Szabo  and  I.  Babuska  presented  at  the  19th 
National  Symposium  on  Fracture  Mechanics,  San  Antonio,  Texas,  June  30- 
July  2,  1986. 

"A  Global  Algorithm  for  Finding  all  Roots  of  an  Analytic  Function  in  a 
Finite  Domain,"  X-R.  Ying  and  I.N.  Katz,  presented  at  the  SIAM  1986 
National  Meeting,  Boston,  Massachusetts ,  July  21-25,  1986. 

"A  Reliable  Argument  Principle  Algorithm  to  Find  the  Number  of  Zeros  of 
an  Analytic  Function  in  a  Bounded  Domain,"  X-R.  Ying  and  I.N.  Katz, 
presented  at  the  Symposium  on  the  Impact  of  Mathematical  Analysis  on 
the  Solution  of  Engineering  Problems,  University  of  Maryland,  September 
17-19,  1986. 


"On  Stress  Analysis  with  Large  Length  Rations,"  by  B . A .  Szabo, 
presented  at  First  World  Congress  on  Computational  Mechanics,  Austin, 
Texas,  September  24,  1986. 

"A  Reliable  Argument  Principle  Algorithm  to  Find  the  Number  of  Zeros  o 
an  Analytic  Function  in  a  Bounded  Domain,"  X-R.  Ying  and  I.N.  Katz, 
presented  at  the  First  International  Conference  on  Industrial  and 
Applied  Mathematics,  June  29-July  3,  1987,  Paris. 
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8 . 3  Visit  and  Seminars  at  Government  Laboratories 


"Advanced  Stress  Analysis  Technology",  by  B.A.  Szabo  and  I.N.  Katz, 
presented  on  September  8,  1977  at  the  Air  Force  Flight  Dynamic 
Laboratory,  Wright-Patterson  Air  Force  Base. 


abstract 


With  one  exception,  all  finite  element  software  systems  have  element 
libraries  in  which  the  approximation  properties  of  elements  are  frozen.  The 
user  controls  only  the  number  and  distribution  of  finite  elements.  The 
exception  is  an  experimental  software  system,  developed  at  Washington 
University.  This  system,  called  COMET-X,  employs  conforming  elements  based 
on  complete  polynomials  of  arbitrary  order.  The  elements  are  hierarchic, 
i.e.  the  stiffness  matrix  of  each  element  is  embedded  in  the  stiffness 
matrices  of  all  higher  order  elements  of  the  same  kind.  The  user  controls 
not  only  the  number  and  distribution  of  finite  elements  but  their  approxima¬ 
tion  properties  as  well.  Thus  convergence  can  be  achieved  on  fixed  mesh. 
This  provides  for  very  efficient  and  highly  accurate  approximation  and  a  new 
method  for  computing  stress  intensity  factors  in  linear  elastic  fracture 
mechanics.  The  theoretical  developments  are  outlined,  numerical  examples 
are  given  and  the  concept  of  an  advanced  self-adaptive  finite  element 
software  system  is  presented. 


2. 


"The  Constraint  Method  for  Finite  Element  Stress  Analysis,"  by  I.N. 
Katz,  presented  at  the  National  Bureau  of  Standards,  Applied 
Mathematics  Division  on  October  19,  1977. 


abstract 


In  conventional  approaches  to  finite  element  stress  analysis  accuracy 
is  obtained  by  fixing  the  degree  p  of  the  approximating  polynomial  and  by 
allowing  the  maximum  diameter  h  of  elements  in  the  triangulation  to  ap¬ 
proach  zero.  An  alternate  approach  is  to  fix  the  triangulation  and  to 
increase  the  degrees  of  approximating  polynomials  in  those  elements  where 
more  accuracy  is  required.  In  order  to  implement  the  second  approach  effi¬ 
ciently  it  is  necessary  to  have  a  family  of  finite  elements  of  arbitrary 
polynomial  degree  p  with  the  property  that  as  much  information  as  possible 
can  be  retained  from  the  pth  degree  approximation  when  computing  the 
(p+l)st  degree  approximation.  Such  a  HIERARCHIC  family  has  been  formulated 
with  p  >  2  for  problems  in  plane  stress  analysis  and  with  p  >  5  for 
problems  in  plate  bending.  The  family  is  described  and  numerical  examples 
are  presented  which  illustrate  the  efficiency  of  the  new  method. 


"The  p-version  of  the  Finite  Element  Method,"  by  I.N.  Katz  and  B.A. 
Szabo,  to  be  presented  at  Air-Force  Flight  Dynamics  Laboratory  Wright- 
Patterson  Air  Force  Base  on  April  23,  1981  (tentative  date). 


abstract 


The  theoretical  basis  of  the  p-version  of  the  finite  element  method  has 
been  established  only  quite  recently.  Nevertheless,  the  p-version  is  al¬ 
ready  seen  to  be  the  most  promising  approach  for  implementing  adaptivity  in 


practical  computations.  The  main  theorems  establishing  asymptotic  rates  of 
convergence  for  the  p-version,  some  aspects  of  the  algorithmic  structure  of 
p-version  computer  codes,  numerical  experience  and  a  posteriori  error  es¬ 
timation  will  be  discussed  from  the  mathematical  and  engineering  points  of 
view . 


U.  Attended  the  meeting  sponsored  by  AFOSR  on  "The  Impact  of  Large  Scale 

Computing  on  Air  Force  Research  and  Development,"  at  Kirtland  Air  Force 
Base,  Albuquerque,  New  Mexico,  April  4-6,  1984. 


5.  Attended  the  meeting  sponsored  by  AFOSR  on  "The  Impact  of 

Supercomputing  on  Air  Force  Research  and  Development  II,"  at  Eglin  Air 
Force  Base,  Florida,  November  3-5,  1986. 


